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REPORT  SUMMARY 


Automation  Technology,  Inc.  reports  that  PE  tests  and  Card 
tests  proceeded  well  during  the  quarter.   Diagnostic  efforts  will  "be 
concentrated  on  "diagnosing"  conversion  problems  through  a  major 
portion  of  the  next  quarter. 

During  this  quarter  the  ILLIAC  IV  Project  has  received  a 
Burroughs  B65OO/B67OO  Computer  at  the  University  of  Illinois.   This 
machine  has  been  obtained  primarily  to  support  the  programmers  devel- 
oping the  ILLIAC  IV  operating  system.   It  is  also  available  to  appli- 
cations programmers  and  to  prospective  users  of  the  ILLIAC  IV  computer. 

The  software  effort  was  directed  toward  integrating  the 
Operating  System  with  the  ILLIAC  IV  simulator  on  the  B650O,  converting 
the  Cockroach- to -Glypnir  translator  to  the  B65OO,  supporting  the 
assembler  (ASK)  on  the  B65OO,  and  simulating  ILLIAC  IV  special  functions 
and  algorithms  written  in  ASK  for  a  mathematical  subroutine  library. 
The  Software  Reference  Manual  is  nearly  complete. 

Several  major  software  efforts  for  the  PDP-11  system  were 
continued  or  completed.   Hardware  delivery  began  on  the  various  seg- 
ments of  the  ARPA  network  terminal  system.   Negotiations  were  completed 
and  agreement  reached  with  Digital  Equipment  Corporation  on  a  co- 
operative research  arrangement  for  joint  development  of  the  ARPA  network 
port  facility. 

Application  efforts  in  Numerical  Analysis  continue  in  solving 
partial  differential  equations,  matrix  inversion  and  linear  algebraic 
equations,  eigenvalues  and  eigenvectors,  estimation  and  filtering, 
and  identification  of  nonlinear  differential  equations. 


Three  seminars  were  offered  during  the  quarter. 

Project  expenditures  and  commitments  through  March,  1971: 

Burroughs  Corporation  $26,723,000.00 

University  of  Illinois  6,8lU,273-83 
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1.   HARDWARE 


1.1  Off -Line  Diagnostics 

1.1.1  PE  Test 

The  analysis  of  the  PE  logic  and  the  development  of  tests  to 
he  used  to  test  the  PE  off-line  through  use  of  the  PE  Exerciser  (PEX) 
is  progressing  "well.   The  Path  Tests  and  Combinational  Logic  Tests  are 
completed  and  in  use.   The  analysis  necessary  for  development  of  tests 
of  the  Control  Logic  is  nearing  completion  and  some  test  sequences  have 
been  defined.   This  area  is  expected  to  be  complete  by  late  summer. 

1.1.2  Card  Test 

The  development  of  tests  to  isolate  stuck-type  faults  in 
active  devices  is  on  schedule  except  as  noted  below.   All  tests  for  PE 
Cards  are  complete ,  have  been  verified ,  and  are  in  use.   Of  the  CU 
Cards,  the  Type  1  boards  are  scheduled  to  be  completed  by  mid- summer. 
The  Type  1  boards  have  been  divided  into  two  classes,  31  Priority 
Boards  and  kO   Non-Priority  Boards.   Twenty-eight  Priority  Boards  and 
21  Non-Priority  Boards  have  their  tests  generated;  however,  these 
tests  are  not  verified,  nor  are  they  in  use,  since  there  is  no  CU  Card 
Tester  operational. 

1.2  Program  Conversion 

Conversion  of  the  programs  used  to  develop  the  diagnostics 
mentioned  above  was  begun  during  the  quarter.   Because  of  the  impor- 
tance of  conversion  to  be  able  to  carry  on  useful  test  generation, 
some  other  scheduled  work  is  being  postponed.   The  level  of  effort  on 
the  PE  Control  Logic  Tests  is  temporarily  reduced  to  analysis  only 
and  no  board  tests  are  being  generated.   The  majority  of  diagnostic 
efforts  will  be  concentrated  on  "diagnosing"  conversion  problems 
through  a  major  portion  of  the  next  quarter. 
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1.3  ILLIAC  IV  Maintenance 

ATI  maintenance  engineers  are  actively  participating  in  the 
on-line  debugging  of  ILLIAC  IV.   Among  the  areas  of  active  participation 
have  been:   PE  Board  Test,  PEM/MLU  Test,  PE  Test,  CU  Test  and  i/O  Test. 
In  addition,  ATI  has  identified  the  long  lead  spare  parts  and  estimated 
the  required  number  of  each.  Under  authorization  of  the  University  of 
Illinois,  ATI  is  obtaining  competitive  prices,  where  possible. 


l.k     Financial 

At  this  point  in  the  contract,  ATI  is  approximately  ten 
percent  (10$)  below  the  budgeted  cost  estimates. 
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2.      SOFTWARE 


2.1  Operating  System 

2.1.1  Operating  System  I 

The  operating  system  is  at  present  being  integrated  with  the 
ILL1AC  IV  simulator  (SSK)  on  the  B650O,  and  has  reached  the  stage  where 
simple  ILLIAC  TV  programs  have  been  run  under  the  control  of  the  ILLLAC 
IV  operating  system.   Instrumentation  and  error  diagnostics  are  now 
being  placed  in  some  of  the  modules  to  allow  the  performance  of  the 
operating  system  to  be  measured  in  detail.   An  Operating  System  Main- 
tenance Manual  has  been  written. 

2.2  Compilers  and  Translators 

2.2.1  Glypnir 

Glypnir,  Version  II,  has  been  undergoing  further  consolidation. 
i/O  for  ILLIAC  IV  (as  opposed  to  Simulator  i/o)  and  Macro  Facilities  in 
the  form  of  a  pre-processor  have  been  provided.   The  Glypnir  compiler 
has  also  undergone  a  basic  measurement  investigation  with  a  view  to 
speeding  it  up.   The  possibility  of  providing  a  facility  for  separately 
compiling  subroutines  has  been  and  is  still  being  studied.   A  Glypnir 
Compiler  Maintenance  Manual  has  been  written. 


2.2.2  Cockroach 

During  the  quarter  the  Cockroach-to-Glypnir  translator  was 
completed  for  a  subset  of  the  specified  language  on  the  B650O  and  con- 
verted to  the  B65OO.   Subroutines  and  functions  will  be  completed  in 
the  beginning  of  the  second  quarter.  User  documentation  of  the  available 
features  was  also  provided.   The  addition  of  an  hourly  employee  to  the 
staff  hastened  the  completion  and  provided  for  more  extensive  debugging. 
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Cockroach  is  now  available  to  users,  and  a  certain  limited  amount  of  user 
support  also  is  available. 

2.3  Assembler 

The  Assembler  (ASK)  is  now  supported  on  the  B65OO.   It 
compiles  at  about  1200  cards  per  minute.   A  plan  for  increasing  its 
compiling  speed  to  2500  cards  per  minute  has  been  formulated  and  is 
being  implemented. 

2.4  Interactive  Communications  and  Graphics 

2.4.1  Interactive  Communications 

During  the  reporting  period  several  major  software  efforts  for 
the  PDP-11  system  were  continued  or  completed. 

1.  A  high-level  language  compiler,  PPESPOL,  was  completed. 
Version  II  is  running  on  the  B65OO. 

2.  Version  I  of  the  ARPA  Network  Terminal  System,  ANTS,  was 
completed.   Basic  portions  of  the  system  were  checked 
out  and  actual  runs  made  on  the  PDP-11. 

3.  Design  of  the  link-up  between  Paoli  and  the  University 
of  Illinois,  via  the  ARPA  network,  using  two  ANTS  systems 
was  completed  and  initiated. 

k.      Final  specifications  were  agreed  upon  by  Burroughs  for 
the  construction  of  an  interface  between  the  B65OO/B67OO 
system  and  the  ARPA  network  IMP. 

During  this  reporting  period,  hardware  delivery  began  on  the 
various  segments  of  the  ARPA  network  terminal  system.   The  basic  PDP-11 
processor,  the  l6K  words  of  memory,  the  real-time  block,  the  high- 
speed paper  tape,  ASR35  operator's  teletype,  four  24-00  Baud  line  inter- 
faces with  adapters,  one  dataset  control  and  interface  for  a  Bell  103A 
dataset,  plus  two  general-purpose  interfaces  to  be  used  to  connect  the 
Gould  Electrostatic  Plotter  to  the  system  and  to  connect  the  Computek 
storage  scope  system  were  received.   The  remainder  of  the  items  ordered 
should  be  delivered  within  the  next  reporting  period. 
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Two  interfaces  between  the  AREA  network  IMP  and  the  PDP-11 
were  received.   Both  interfaces  were  debugged,  checked  out,  and  in- 
stalled in  the  University  of  Illinois  PDP-11  and  the  Paoli  PDP-11. 
Success  in  getting  on  the  network  now  hinges  on  the  completion  of 
system  software  development. 

Negotiations  were  completed  and  agreement  reached  with 
Digital  Equipment  Corporation  on  a  cooperative  research  arrangement 
for  joint  development  of  the  ARPA  network  port  facility. 

2.4.2  Graphics 

No  work  was  done  in  the  graphics  area  during  this  reporting 
period  because  of  total  lack  of  personnel.   An  effort  was  begun,  toward 
the  end  of  the  period,  to  acquire  a  full-time  professional  in  graphics 
who  is  expected  to  join  the  Center  during  the  next  period. 

2. 5  Mathematical  Subroutines  Special  Function  Library 

This  quarter's  work  has  continued  on  the  ILLIAC  IV  Special 
Functions  and  Algorithm  subroutine  library.   The  following  set  of 
routines  written  in  ASK  have  been  successfully  simulated  by  the  current 
version  of  the  simulator  and  the  results  are  good  for  15  significant 
digits : 

64-bit  mode:   *sine  and  cosine 

•  tangent  and  cotangent 
•natural  logarithm 

•  square  root 

•arctangent  (l  argument) 
•arctangent  (2  arguments) 
•hyperbolic  sine  and  cosine 
•hyperbolic  tangent 
•gamma  range ( 0,1) 

•  gamma  range (0,<») 

32-bit  mode:   *sine  and  cosine 

•tangent  and  cotangent 
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•exponential 
•natural  logarithm 
•square  root 

Algorithmic  library: 

•addition/ subtraction  of  matrices  stored  straight 
or  skewed 

•multiplication  of  matrices  stored  straight 

•transpose  of  a  matrix  stored  straight 

At  present  work  is  being  completed  for: 

64-bit  mode:    'gamma  (X)  0  <  X  <  57 

57  <  X  <  oo  auid 
X  <  0 
•error  function:  erf  (x)  and 

LN  gamma  (X) . 

For  matrix  operations  the  set  is  being  expanded  to  rectangular 
matrices  stored  striaght  and  skewed.   The  same  set  will  also  be  made 
available  in  Glypnir.   These  routines  will  continue  to  be  refined. 
Maintenance,  distribution  and  programming  assistance  in  using  the 
mathematical  routines  is  available. 

2.6  B5500/B6500  Operations 

During  this  quarter  the  ILLIAC  TV   Project  has  taken  delivery 
of  a  Burroughs  B650O  computer  at  the  University  of  Illinois.   This 
machine  has  been  rented  primarily  to  support  the  programmers  developing 
the  ILLIAC  TV   operating  system.   It  is  also  available  to  applications 
programmers  from  the  Center  for  Advanced  Computation  and  future  users  of 
ILLIAC  IV.   The  machine  was  installed  in  the  basement  of  the  Coordinated 
Science  Laboratory  and  became  operational  around  March  1. 

2.7  Software  Documentation 

The  Software  Reference  Manual  is  now  at  900  pages,  largely 
completed.   Its  one  deficiency  is  an  Assembler  users  manual,  which  is 
at  present  about  half  completed. 


3-   APPLICATIONS 

3.1  Numerical  Analysis 

3.1.1  Partial  Differential  Equations 

3.1.1.1  Numerical  Solution  of  Problems  in  Hydrodynamics 

The  results  of  numerical  experiments  with  the  Brailovskaya, 
Dufort-Frankel,  Cheng-Allen,  Crank-Nicholson  and  Lax-Wendroff  finite 
difference  schemes  on  the  Burgess  equation  have  been  presented  in  a 
formal  report  which  is  now  in  the  process  of  being  printed  [  1] . 

A  code  is  being  written  and  debugged  for  the  two  dimensional 
subsonic  and  transonic  flow  around  a  circular  cylinder,  as  the  next 
step  in  exploiting  fully  the  techniques  of  parallel  processing  in 
the  numerical  computations  of  three-dimensional  fluid  and  gas  flow 
problems.   In  a  polar  co-ordinate  system  the  governing  equations  in 
non-dimensional  form  are: 


Continuity    r—   +■  py  .  V  =  0: 
Dt 


Momentum     p  ~   +.  yp  =  <7   ~  (4/3  W-V-yxyxV); 


DS    r~  /      \  Moo       s/7  Moo   2  /  /  N 
Energy       p  gg  =  Jy    (7  -  1)  g  <£>  i  7  ^^  y  (p/p) 


V  =  velocity,  p  =  pressure,  S  =  entropy,  7  =  ratio  of  specific  heats, 
M»  =  free  stream  Mach  number,  $  =  dissipations  terms,  Re  =  Reynolds 
number,  Pr  =  Prandtl  number. 

The  computational  region  comprising  the  flow  field  is  divided 
into  a  suitable  number  of  mesh  points,  which  is  a  function  of  computer 
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storage  capacity,  and  the  dependent  variables,  namely  the  pressure  p, 
radial  and  tangential  components  of  the  velocity,  i.e.  u,  v,  and  tem- 
perature, T,  are  calculated  at  each  mesh  point.   The  equation  of  state 
for  a  perfect  gas  is  assumed.   Initially  the  gas  is  assumed  to  be 
flowing  with  a  uniform  velocity  and  the  physical  properties  are  assumed 
uniform  at  each  mesh  point. 

To  gain  flexibility  and  save  time  on  runs  while  the  program 
is  being  debugged,  the  viscous  terms  have  been  suppressed  in  the 
calculation,  so  that  the  flow  is  treated  as  being  inviscid.   At  time 
t  =  0,  the  condition  of  zero  radial  velocity  is  imposed  on  the  cylinder 
boundary,  a  suitable  time  step  is  selected,  and  the  finite  difference 
representation  of  the  governing  equations  is  used  to  calculate  the 
values  of  the  dependent  variables  at  succeeding  times.   The  finite 
difference  scheme,  chosen  for  these  experiments  is  the  two-step, 
second  order  accurate,  Richtmyer  variation  of  the  Lax-Wendroff  scheme 
[  2] ,  which  has  been  previously  tested  with  the  Burgess  equation  [  1] . 

3.1.1.2  Algorithm  Development 

During  this  period  several  numerical  methods  have  been  in- 
vestigated in  cooperation  with  AMDCO  Oil  Company  to  solve  the  elliptic 
type  partial  differential  equations.   The  semi-iterative  block  Jacobi 
method  and  ADI  method  were  adopted  and  both  were  implemented  in  Glypnir. 
The  first  method  worked  correctly  and  was  made  into  a  general  type 
subroutine.   An  iterative  step  takes  about  600  microseconds  on  ILLIAC  IV 
in  the  case  of  21  x  21  mesh  points.   The  second  method  is  still  being 
implemented.   The  storage  scheme  for  this  method  is  more  complicated 
than  the  first  method. 

3.1.2  Matrix  Inversion  and  Solution  of  Linear  Algebraic  Equations 

During  this  quarter,  a  program  to  solve  linear  equations  and 
perform  matrix  inversion  has  been  written  in  ASK,  and  is  presently 
being  debugged.   It  uses  the  Gaussian  Elimination  method  and  is  designed 
to  handle  full  matrices  up  to  630  x  63O  for  Gaussian  elimination  and 
up  to  ^75  x  ^75  for  inversion.   These  large  matrices  are  not  core 
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containable  and  are  therefore  "broken  down  into  blocks  of  6k   columns, 
each  block  being  processed  independently  of  the  others  as  far  as 
possible. 

The  algorithm  works  in  two  stages,  first  the  matrix  A  is 
decomposed  into  the  product  of  lower  and  upper  triangular  matrices, 
L  U. 

The  second  stage  uses  this  decomposition  to  operate  on  a 
right  hand  side  to  solve  Ax  =  b,  on  a  unit  matrix  to  form  A   or  on 
an  arbitrary  matrix  B  to  form  A  B.   To  reduce  the  round-off  error, 
'partial  pivoting'  is  employed  (see  e.g.  Wilkinsen,  J.  H. ,  The 
Algebraic  Eigenvalue  Problem) ,  also  the  inner  products  of  vectors  are 
accumulated  using  double  precision  arithmetic. 

Work  is  just  beginning  on  the  Conjugate  Direction  method 
which  will  be  implemented  in  Glypnir.   This  should  prove  an  efficient 
algorithm  for  solving  Ax  =  b  when  the  matrix  A  is  not  given  explicitly 
but  a  subroutine  which  computes  Ax  given  the  vector  x  is  available. 

3.1.3  Eigenvalues 

Jacobi's  Method  for  finding  eigenvalues  and  eigenvectors  of 
real  symmetric  matrices  (including  complex  Hermitian  matrices:  Let 

A  =  B  +-  iC  be  a  complex  Hermitian  matrix  where  B  is  real  symmetric 

t  t 

(B  =  B  )  and  C  skew- symmetric  (C  =  -C  ),  then  the  real  symmetric 

2n  x  2n  matrix 


a1- 


B 


-C 


can  be  an  input  to  the  Jacobi  algorithm)  has  been  coded  in  ASK  and 
satisfactorily  tested  on  the  B5500  ILLIAC  IV  simulator. 

Eberlein's  Method  (jacobi-like  algorithm)  for  normalizing 
real  non-symmetric  matrices  has  also  been  satisfactorily  tested. 
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The  results  of  both  algorithms  agree  with  the  literature  up  to 
Ik   significant  digits.  An  error  of  order  0  (10   )  is  due  to  conversion 
routines  from  internal  machine  representation  of  floating  points  to 
their  external  decimal  representation. 

To  have  a  means  of  comparison,  two  well  established  algorithms 
by  H.  Rutishauser  [3]  and  P.  J.  Eberlein  [h]   have  been  combined  in  one 
B67OO  ALGOL  program.   It  is  in  a  running  condition  and  exists  under  the 
file  name  (completed  version)  WINFRIED/BEKMARD/PROC.   This  file  is  a 
separately  compiled  procedure.   The  name  of  the  procedure  is  JACEIG 
(N,A,TS,D,R0T,TMX). 

In  using  this  procedure  the  user  has  an  option 

1)  He  may  write  his  own  main-program  for  i/O  and  procedure 
call,  but  he  then  has  to  go  through  the  B67OO  BINDER  for 
interrelation  between  his  file  and  the  above  file  or, 

2)  He  can  use  the  file  BERNHARD/jACOBI/EIGEN  which  does  the 
i/O  and  calculation  for  him. 

For  case  (2)  the  following  control  cards  and  input -data  are 
required : 

5  EXECUTE  BERNHARD/jACOBI/EIGEN 

5  BCL  EBOR 

<integer  1>,   <integer  2>, 

matrix-elements   in  real. 

8  END 

where  5  is  an  invalid  character  like  multipunch  1,2,3 

<integer  1>  =  order  of  matrix 
<integer  2>  =  TMX  with 

TMX  =  0  no  eigenvectors  are  calculated  and  printed 
TMX  >  0  the  right  eigenvectors  are  calculated 

and  printed 
TMX  <  0  the  left  eigenvectors  are  calculated 
and  printed. 
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The  matrix  for  which  the  eigenproblem  is  to  be  solved  can  be 
either  real  symmetric  or  real  non-symmetric.   If  real  symmetric  then 
only  the  eigenvalues  are  printed  in  decreasing  order.   If  the  eigen- 
vectors are  to  be  printed  (TMX  #  0),  then  this  will  leave  the  same 

order  as  the  eigenvalues,  i.e.,  to  \.    corresponds  TS .  in  the  order  of 

J  J 

their  appearance. 

If  real  non-symmetric  the  total  eigenvalue  matrix  is  printed. 
If  the  eigenvalues  are  complex  then  the  a. .  occupy  the  diagonal  while 
the  imaginary  part  occupies  the  off  diagonal  portions . 

If  a. .  +  i  a. .  is  the  eigenvalue,  the  corresponding  eigen- 
11  —    ij 

vector  t.  +  i/+  \    appears  in  column  (row)  i  and  column  (row)  j. 
J 

3.1.3.1  Eigenvalues  and  Eigenvectors  of 
Symmetric  Tridiagonal  Matrices 


3.1.3.1.1  QR  Algorithm  with  Origin  Shifts 

The  QR  algorithm  for  finding  eigenvalues  of  symmetric  tri- 
diagonal matrices  has  been  written  in  XALGOL ,  debugged  and  tested  for 
different-size  matrices.   A  short  description  of  the  algorithm  follows. 


A  is  a  symmetric  tridiagonal  matrix  with  diagonal  elements 
diagonal  elements  b.  .   A] 
A  ,  (t  =  1,  2,  ...),  may  be  expressed  as 


C.  and  subdiagonal  elements  b. .   Any  symmetric  tridiagonal  matrix 


\  ■  «tR, 


(i) 


where  Q     is   orthogonal   and  R     is   upper  triangular  of  the   form 


pl  ql  rl 

p2  <i2  r2 


\       \. 


V\  *V2 

5n-l 


n 
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Let  At+1  ■  <*t  At  \     then     Vl  ■  Rt  V 


(2) 


For  large  t,  A  approaches  a  diagonal  matrix,  Diag  [X.],  where 
X. 's  are  the  eigenvalues  of  the  matrix,  A. 

To  accelerate  convergence  of  this  algorithm,  we  let 


A'   =  A     -KI 

X>  "G  "G 


where  K     is  the  eigenvalue  of  the  2x2  matrix 


.00 

'n-1 


n-1 


which  is   closer  to   C 


(t) 


,(t) 

n-1 


n 


n 


Then  we   decompose   A'    =  QL  R.    and  find 


W  =  V»t  +  V 


By  repeating  these  two  steps,  the  eigenvalues  of  the  matrix,  A,  are 
immediately  available. 


A  document  describing  the  algorithm  with  flow  chart  is  "being 


written. 


3.1.3.1.2   Inverse  Iteration  to  Find  Eigenvectors 

A  program  to  calculate  the  eigenvectors  of  symmetric  tri- 
diagonal  matrices  with  given  eigenvalues  has  been  written  in  XALGOL, 
debugged  and  tested.   The  algorithm  is  as  follows: 

A  is  a  symmetric  tridiagonal  matrix  with  c.  as  its  diagonal 
elements  and  b.  as  its  subdiagonal  elements.   Its  eigenvalues  are  given 
rather  accurately  and  arranged  in  descending  order. 

Let  vector  V  be  taken  as  the  initial  vector.   We  solve 
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(A  -  \I)  X  =  V  (1) 

followed  by 

(A  -  \I)  Y  =  X  (2) 

then  Y  is  the  eigenvector  corresponding  to  eigenvalue  X. 
To  make  the  algorithm  more  effective,  we  let 

(A  -  XI)    =  LU  (3) 

where  L  is  unit  lower  triangular  and  U  is  upper  triangular,  set  V  =  Le 
where  e  =  (l,  1,  ...,  l),  then 

U  X  =  e .  ( k ) 

Hence  X  can  he  determined  by  back-substitution.   The  matrix  U  is 
determined  by  Gaussian  elimination  with  partial  pivoting,  i.e.  eliminate 
the  variables  in  their  natural  order,  but  select  the  row  having  the 
maximal  coefficient  of  X.  as  the  pivotal  row  in  the  i-th  stage. 
Similarly,  we  can  find  the  vector  Y. 

A  document  describing  the  algorithm  with  flow  chart  and 
numerical  examples  is  being  written. 

3.1.U  Polynomial  Root  Finder 

Actual  coding  of  a  root  finder  has  not  been  done  in  ASK. 
However  ALGOL  codings  have  been  debugged  on  the  B5500.   Presently 
changes  are  being  made  in  these  programs  to  run  on  the  B6500. 

In  general,  the  algorithm  has  been  outlined  to  be: 
(l)  evaluate  polynomial  at  various  intervals  to  determine  whether  a 
change  of  sign  takes  place,  (2)  shorten  these  intervals  as  much  as 
possible,  (3)  proceed  on  new  intervals  with  Newton-Rapheson  iterations. 
Also  an  ALGOL  program  using  Sturm  sequences  to  evaluate  intervals  is 
presently  being  debugged  on  the  B65OO . 


-15- 


3.1.5  Identification  of  Non- linear  Differential  Equations 

Much  is  known  about  the  numerical  solutions  of  differential 
equations,  however,  almost  nothing  is  known  about  the  reverse  problem, 
i.e.,  given  some  observed  function  of  time,  x(t),  can  we  find  a 
differential  equation  of  which  x(t)  is  a  solution? 

Such  problems  exist  in  economics,  chemistry,  medicine  and 
many  other  fields,  the  equations  involved  often  being  nonlinear.   In 
the  past  solution  has  been  restricted  to  linear  problems  or  problems 
involving  the  estimation  of  only  a  few  parameters.   This  has  been  caused 
not  so  much  by  lack  of  solution  techniques,  but  by  inadequate  computing 
power.   ILLIAC  IV  can  alleviate  this  problem  through  its  high  computing 
speed  and  the  fact  that  the  problem  can  utilize  the  array  processor 
efficiently. 

Since  many  different  differential  equations  may  have  the  same 
solution,  it  is  necessary  to  restrict  the  problem  to  invariant 
equations  of  the  form: 

F(x,  x',  x",  ...,  x(m),  fx,  ...,  fn)  =  0  (1) 

where  F  is  assumed  to  be  known, 

x(j)  =  J^     x(t)  and 
dtJ 

f.  =   f.(x^i')  are  unknown  functions 
i    i 

c 

or  x  or  one  of  its  derivatives.   The  problem  is  to  develop  an  algorithm 

to  determine  the  f.  . 

l 

E.G.     x"  +  f(x)  =  0 

If  x  were  given  as  x(t)  =  sin  t,  we  could  derive  f(x)  =  x.   Let  x(t)  be 
the  observed  function,  defined  on  the  interval  0  <  t  <  T  . 

The  solution  to  (l)  depends  on  the  functions  f .  and  the 
initial  conditions  of  x  and  its  first  n  -  1  derivatives  at  t  =  0. 
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I.e.,  x  -  x   (.x_,    x_,    x_,    . ..,    Xq  ,    f  ,    • ..,    f   ) 


or  more  briefly  x  =  x(x„,    f) 


We  define  the  error  functional 

E  Qv  £)  =  II  *  -  x  (*0>  £)  II2 


where 


II  y  II2  =  Mt)  [y(t)]2  dt 
J  o 

for  some  weighting  function,  w.   E  is  a  measure  of  the  difference  between 
observed  function  and  the  computed  one  for  some  x_  and  f . 

The  problem  thus  reduces  to  finding  x~  and  f  such  that  E  (x_,  f) 
is  minimized. 

For  the  numerical  solution,  x(t)  is  assumed  to  be  observed  at 
discrete,  equal  intervals  and  the  integrals  replaced  by  summations.   It 
is  also  assumed  that  these  observations  are  contaminated  by  random 
"noise"  so  that  minimizing  E  constitutes  a  least  squares  solution  to  the 
problem. 

The  functions  f .  are  approximated  over  the  interval  on  which 
they  are  to  be  identified,  [a,b],  by  a  linear  combination  of  fixed 
orthonormal  functions 

n 
f(x)  m   Z   q.  $.  (x) 

where  <  $.,  <$>.  >_  =  8.  .,  the  Kroneuker  delta  function,  the  inner  product 
I'  3     1    10  . 

<  g,  h  >  being  given  by  <  g,  h  >  =  /  g  (x)  h  (x)  dx.   In  this  imple- 

mentation  the  functions,  $>.,  will  be  splines. 

J 

The  Solution  Algorithm 

For  simplicity,  we  consider  the  problem 
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n 
Let  f(x)  =   Z   q.  $.  (x)  and  x  (0)  =  x  =  q^.   Define  the  solution  of 
i=l 

x'    =  f(x) 
x  (0)  =  xQ 

to  be  x  (t).   The  dependence  on  x  being  implicit. 
Then  E(f)  =  E(q)  =  |  |?x|  |2 

f   T   -  2 

=  /   [x(t)  -  x  (t)T  dt. 

J   Q 

When  E(q)  is  minimized,  V  q  E  =  0 

where  (V  q  E) .  =  ~-  E.    i  =  0,  1,  ...,   n. 

A  vector  q  such  that  V  q  E  =  0  may  be  found  iteratively  using 
the  Conjugate  Gradient  Method  [1]  : 


Ti+1    Ti    n  n 
where:    p    =  -  (V  q  E)    , 

P  ,  =  -  (V  q  E)(n+1)  +  b  p   ,  (2) 

n+1  '  n  n 

(V  q  E)(n+1)  =  V  q  E^), 

<(VqE)(»11,(A«)MP>? 
b  ' 

(Ae)m=AmVi) 

and  v  q  E  is  the  matrix  of  second  derivatives: 
(^  q  E) 


Id  '"  ^  c*L- 
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The  superscripts  indicate  the  iteration  count.   The  inner  product 
<  a,b  >  "being  defined  by 

n 
<  a.b  >  =   E    a.  "b.  . 
2    i=0    x     x 

At  each  iteration,  c  is  chosen  so  as  to  minimize  E(g  +  cp  )  over  all  c. 

n  n     n 

The  calculation  of  -sr —  E  . 

^-E  =  ^  jw(t)  [x(t)  -xf(t)]2dt 

f   T  ^f  ft) 

=  -2  Jw(t)  [x(t)  -  xf(t)]  .   ^  [t)      dt  (3) 

oX* 

if  we  abbreviate  ^ —  by  6  x.,  then  5  x.  satisfies  the  so-called 

oq..       i         i 

i 

"sensitivity  equations" 

dt  5  xo  =  ^  (xf}  '  5  V  6  xo  (0)  =  1 

and      —  S  x.  =t-  (xj  .   5  x.  +  0.  (xj,  5  x.  =0,  i  =  1,  2,    ...,    n 
dt     i   dx   f       i    if      l 

representing  the  sensitivity  of  the  equation  x'  =  f(x)  to  changes  in  xn 
and  f  respectively. 

If  we  define  $  (x)  =  0,  these  may  be  combined  to: 

dt  5  xi  =  f  (xP  •  8  xi  +  *i  (xP  1  -  0,  1,  ....  n  .    (U) 


To   determine  v     q  E  we  note  that, 


2  2  -  T 

■'w(t)    [n(t)    -  n-  (t)]2     dt 


M±   oXL.         aq±    o^       J0  f 

=  • 2  J{(t)  ^  xf(t)  ^  xf(t)   dt 

,-  T  2 

-   2      /w(t)    [n   (t)    -  xf(t)]    x   ^/^    (xf)    (t)    dt 
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=  2      /w(t)    5  n.    6  n     dt 
J  n  x  J 


-  2     /w(t)    [n(t)   -  xJt)]    6  nf .   dt  (5) 

2                d2 
where  we  define  5  n.  .  =  « e —  (x„)    (t)    . 

o 

By  differentiating  (t)  with  respect  to  q.;  5  n. .  may  be  shown  to 
satisfy  the  equation: 

d  *2      ^f  »    » 
-rr   o  n.  .  =  — =■  5  n.  on. 
dt    ij    ^2     i    j 

df   o2 

+  t—  on.. 

dn       1J 

+  4-  «.  S  n.  +  |-  <D.  &  .  (6) 

dn   i    j    dn   J  ni 

Since  we  require  only  v  q  E  explicitly,  we  need  calculate  only 

n    2 
the  components  ¥.  =   Z   5  n. .  p.  . 

1    j=0     1J   J 
Multiplying  (6)  by  p.  and  summing  over  j  gives: 

J 

2 

*  y.  =1|  5  n.  (  £  6  n.  p.) 
dt   i   ^2     i         3   V 

dn   l 


+  k  \  (Z  5n3pP 


+  <$;  z  *j  V 8  ni  (7) 

where  all  summations  are  from  j  =  0  to  n.   This  reduces  the  number  of 
differential  equations  to  be  solved  from  (n  +  l)  to  n  +  1  . 
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Using  (5)  we  obtain 

[^  q  E  p].  =2      /w(t)   S  n.  x  (  I  6  n.  p.)  dt 

-  2   /w(t)  [n  (t)  -  x  (t)]  1   (t)  dt  . 

J   Q 

Preliminary  experiments  indicate  that  this  algorithm  is  rapidly 
convergent,  only  5-10  iterations  "being  required  when  a  good  estimate  of 
f  is  available. 

Implementation  on  ILLIAC  IV 

A  specialized  compiler  is  needed  to  perform  the  following 
tasks: 

1.  Accept  the  form  of  the  equation  to  be  identified. 

2.  By  using  symbolic  differential  techniques,  form  the 
sensitivity  equations. 

3-1.6  Estimation  and  Filtering 

3.1.6.1  Non-linear  Least  Squares  Problem 

Work  is  proceeding  on  the  parameter  estimation  of  non-linear 
least  squares  problems.   The  numerical  algorithms  to  be  investigated  are: 

1.  modified  Gauss  method 

2.  gradient  method 

3.  variable  metric  method 
k.  factorization  method 

The  implementation  of  algorithm  (l)  has  been  successfully 
completed  and  the  ASK  code  for  (2)  is  proceeding  for  which  more  accurate 
results  are  expected.  Work  has  also  begun  on  the  variable  metric  method. 

3-1.6.2  Numerical  Solution  of  the  Non- linear  Matrix  Riccati  Equation 

An  eigenvector  solution  of  the  matrix  Riccati  equation  relating 
to  optimal  control  theory  has  been  studied  during  this  quarter.   For  the 
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cases  in  which  the  matrix  of  coefficients  resulting  in  linear  dynamic 
systems  can  be  transformed  to  a  diagonal  matrix,  there  is  a  straight 
forward  method  for  constructing  a  similarity  transformation  whether  or 
not  the  eigenvalues  are  distinct. 

3.2  Linear  Programming 

This  quarter  has  seen  the  switch  to  the  B65OO  computer  instal- 
lation.  Regrettably,  the  transition  involved  a  degradation  of  the  B5500 
assembler- simulator  support  for  a  prolonged  period,  without  comparable 
facilities  being  available  elsewhere.  This  seven  week  stretch  was  used 
to  bring  LP-system  documentation  into  usable  shape. 

The  simulation  that  could  be  performed  was  directed  wholly 
toward  the  INVERT  package.   There  are  indications  that  simulation  will 
not  permit  sufficient  code  "execution"  to  debug  the  inversion  routines, 
but  work  is  continuing. 

B5500  SETUP  routines  have  been  run  on  the  new  computer,  but 
require  modifications  and  extensive  updating  to  reflect  the  increasing 
load  of  SETUP  that  the  ILLIAC  IV  will  perform.   Experience  with  our  B65OO 
indicates  that  only  the  most  minimal  processing  should  remain  on  this 
easily  overtaxed  facility.   The  future  will  see  consideration  being  given 
to  transferring  more  functions,  including  those--such  as  sorting- -which 
are  ungainly  and  difficult  to  code  for  the  ILLIAC. 

The  problem  of  removing  the  solution  and  associated  information 
is  now  being  resolved.   Initially  these  SETDOWN  procedures  will  be 
executed  as  a  separate  program  as  it  is  estimated  that  abnormal  termina- 
tions will  be  frequent.   Again,  as  with  SETUP,  the  majority  of  the  work — 
such  as  rescaling  and  transforming  data  to  B-65OO  representation — will  be 
executed  on  ILLIAC  IV. 

Documentation  was  reviewed,  rewritten,  and  reorganized  with 
the  intent  of  improving  the  clarity  of  the  system's  interconnections  and 
interpre sumptions.   Not  surprisingly,  discrepancies  were  found.   Further, 
refinement  will  be  necessary,  but  will  be  deferred  until  after  the  code 
has  been  run  on  the  prototype  ILLIAC. 
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3. 3  Long  Codes 

An  analysis  of  the  noise  vectors  used  for  previous  studies 
disclosed  an  alarming  degree  of  correlation  between  the  supposedly- 
independent  random  vectors.   Consequently,  an  improved  pseudo-random 
number  generator  of  proven  period  and  potency  was  substituted  for  the 
previously  used  ad  hoc  generator.   All  identification  algorithms  tested 
were  more  effective  when  used  with  the  less  correlated  noise  vectors. 

A  Kalman  filter  algorithm  is  being  implemented  to  allow  testing 
of  the  parameters  obtained  from  the  identification  algorithms  without 
knowledge  of  the  true  values  of  the  parameters.   If  the  innovation 
sequence  of  a  Kalman  filter  using  the  estimated  parameters  as  a  model  and 
using  the  observed  values  as  input  is  a  white  noise  sequence,  identifica- 
tion can  be  considered  successful. 

The  mest  recently  developed  correlation  approach  for  identifi- 
cation showed  considerable  promise  before  its  use  was  suspended  for 
recoding  to  run  on  the  B65OO.   The  equations  for  this  scheme,  like  most 
other  relations  that  occur  in  identification  processes,  include  both  the 
covariance  matrix  of  the  plant  noise  and  the  covariance  matrix  of  the 
observation  noise,  which  are  usually  unknown.   Manipulation  of  the 
correlation  algorithm  expressions  indicates  that  it  will  be  possible  to 
obtain  not  only  estimates  of  the  parameters  of  the  system,  but  also 


estimates  of  the  two  covariance  matrices.   It  may  be  necessary  to  require 

2 
the  correlation  matrix  of  the  observation  noise  to  be  in  the  form  a  I 

(I  =  Identity  matrix),  but  this  restriction  can  probably  be  relaxed. 


3.^  Signal  Processing 

A  Glypnir  subroutine  has  been  written  which  computes  the  fast 
Fourier  transform  of  any  real  vector  of  length  N,  where  N  is  a  power  of 
two  and  N  >  6k. 

A  previously  written  Glypnir  autocorrelation  program  has  been 
put  into  the  form  of  a  subroutine. 

Programs  are  being  written  for  use  on  the  B65OO  computer  which 
will  be  used  to  thoroughly  test  previously  written  ALGOL  procedures  and 
Glypnir  subroutines. 
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Work  is  continuing  with  Amoco  Production  Company  on  implementing 
some  signal  processing  and  partial  differential  equations  algorithms  on 
the  ILLIAC  IV. 

3-5  Education 

3.5.1  ILLIAC  IV  Seminars 

In  addition  to  the  graduate  course  on  ILLIAC  IV,  three  seminars 
were  offered  this  past  quarter:   A  one-day  seminar  was  given  on  January  5> 
1971-   The  outline  follows: 

Seminar  on  ILLIAC  IV 

Background  9:00  -  9: 30  am 

Buffer,  Pipeline  and  Multiprocessor 

Hardware  Structure  9^30  -  12:00  n 

A  General  Description  of  Array 

A  Sample  Problem 

A  More  Detailed  Description  of  Array 

A  General  Description  of  I/O  System 

Test/Repair  Equipment  and  Diagnostics  1:30  -  2:00  pm 

Physical  Characteristics 

Software  2:00  -  h:  30  pm 

Programming  Languages  --  ASK,  GLYFNIR,  FORTRAN 

Operating  System 

Utilities 

Applications  U:  30  -  5:00  pm 

On  January  11-15,  a  one-week  seminar  was  given  to  Pan  American 
Petroleum.   On  March  15-19;  another  one-week  seminar  was  given  to  employees 
at  Ames  Research  Center,  Moffett  Field,  California,  where  ILLIAC  IV  is  to 
reside.   Additionally  a  one  and  one-half  hour  overview  on  ILLIAC  IV  was 
presented  to  about  300  Ames  employees.   The  outline  for  the  one-week 
seminar  follows. 
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A  Seminar  on  the  ILLIAC  IV  System 
I.   Background  --  Conventional  and  Unconventional  Organizations 

A.  Why  ILLIAC  IV? 

1.   Conventional  Organization 

B.  How  to  speed  up  the  Operation-Design  Philosophies 

1.  Overlap 

a.   Buffer 
"b.   Pipeline 

2.  Replication 

a.   General  Multiprocessor  --  Distribute,  Memory,  ALU,  CU 
i.   Recentralize  Memory 
ii.   Recentralize  ALU 
iii.   Recentralize  CU  --  basis  for  ILLIAC  IV 

3.  Both 

C.  ILLIAC  IV  is  a  Vector  Processor 

II.   Hardware  Structure 

A.  Organization  Chart 

B.  ILLIAC  IV  Array  --  General  Description 

1.  Control  Unit  (CU) 

2.  Processing  Element  (PE) 
3-   Data  Paths 

a.  Control  Unit  Bus  (CU  Bus) 

b.  Common  Data  Bus  (CDB) 

c.  Routing  Network 

d.  Mode  Bit  Line 

C.  An  Illustrative  Problem 

1.  DO  10  I  =  1,  N 

10  A(l)  =  B(I)  +  C(I) 

a.  N  =  6k 

b.  N  <  6k 

c.  N  >  6k 

2.  DO  10  I  =  2,  6k 

10  A(I)  =  B(l)  +  C(l-l) 

a.  Skew  at  Compile  Time 

b.  Skew  at  Execution  Time 
3-   DO  10  I  =  2,    6k 

A(l)  =  B(l)  +  A(l-l) 

D.  ILLIAC  IV  Array  —  A  More  Refined  Description 

1.  PE 

a.   RGD 

2.  PU 
3-   CU 

a.  ADVAST 

b .  FINST 

c.  MSU 

d.  TMU 

e.  ILA 

E.  Another  Illustrative  Problem 

1.   Laplace's  Partial  Differential  Equation  describing  the 
Steady-State  heat  distribution  on  a  slab 

F.  Data  Allocation 
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ILLIAC  IV  I/O  System 

1.   I/O 

Subsystem 

a. 

CDC 

b. 

BIOM 

c. 

IOS 

2.   Disk  File  System 

3-  B6500 

a. 

CPU 

b. 

Memory 

c. 

Multiplexor 

d. 

Peripherals 

i.   Remote  Terminals 

e. 

Laser  Memory 

f. 

ARPA  Network 

III.   Configurations  at  CAC  and  Paoli 
IV.   Diagnostics  and  Test/ Rep air  Equipment 

A.  IDIAP 

B.  PEX 

C.  PEMX 

D.  CUCT 

E.  Some  Physical  Characteristics 

1.   Slides 

V.   Programming  Languages 

A.  ASK 

1.  Background,  Review,  Notation,  Conventions 

2.  Sample  Problems 

a.  Summing  an  Array  of  Numbers 

b.  Finding  the  Maximum  Value  in  an  Array  of  Numbers 

c.  Matrix  Multiplication 

i.   Skewed  Storage 

d.  Temperature  Distribution  on  a  Slab 

i.   Case  1:   one  temperature  value  per  PEM 
ii.   Case  2:   eight  temperature  values  per  PEM 

B.  GLYPNIR 

C .  FORTRAN 

VI.   Operating  System 

A.  ICL 

B.  Utilities 

VII.   Some  Applications 
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k.      ADMINISTRATION 

k.l     Administration  and  Services 

4.1.1  Financial  Report 

Actual  expenditures  and  obligations  for  January-March  1971: 

January       February      March 

Burroughs  Corporation     $312,000.00   ($109,000.00)     N.A. 
University  of  Illinois     183,220. 3^     21+0,118.1+2    $168,093.51 

Expenditures  and  obligations  to  date  through  March  1971 : 

Burroughs        $26,723,000.00 
University        6,8ll+,273.83 

Budgeted  expenditures  -  3rd  quarter,  fiscal  1971* 

January       February      March 

Burroughs  $567,900.00    $298,6^3-00   $252,758.00 

University  215,302.00     215,302.00    215,302.00 


Monthly  status  report  not  received  from  Burroughs'  as  of  h/^O/'Jl, 
therefore  figures  for  March  are  not  available. 
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